Collisional relaxation of two-dimensional gravitational systems 
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Abstract 

Systems with long range interactions present generically the formation of quasi-stationary long-lived 
non-equilibrium states. These non-Boltzmann states relax to Boltzmann equilibrium following a dy- 
namic which is not well understood. In this letter we derive a simple approximate kinetic equation 
for the relaxation process in a two-dimensional inhomogeneous self-gravitating particle system, ob- 
taining a Fokker - Planck equation for the velocity distribution with explicit analytical diffusion 
coefficients. Performing molecular dynamics simulations using the full dynamics and comparing 
them with the evolution predicted by the Fokker - Planck equation, we observe a good agreement 
with the model for all the duration of the relaxation, from the formation of the Quasi-Stationary 
state to thermal equilibrium. During all this process we observe a scaling of the relaxation time 
proportional to the number of particles in the system. 

PACS numbers: 04.40.-b, 05.70.Ln, 05.70.-a 



Systems of particles with long range interactions are 
those which inter-particle potential at large separation 
decays slower than the dimension d of space, i.e., v(r — ¥ 
oo) ~ 1/r 7 with 7 < d. There are many examples in 
nature such as self-gravitating systems in the cosmologi- 
cal and astrophysical context (the large structure of the 
universe, galaxies, etc), interaction between vortices in 
two-dimensional hydrodynamics, cold classical atoms or 
capillary interactions between colloids or granular media 
(for a review see e.g. These kinds of systems present 
very particular properties in thermal equilibrium, such 
that negative micro-canonical specific heat or inequiva- 
lence of statistical ensembles. Their dynamics is also pe- 
culiar compared to short range systems: in a first stage 
there is the generic formation in a few characteristic times 
Tdyn of a long-lived non-equilibrium state — during the 
so-called violent relaxation process. A typical example of 
such quasi-stationary states (hereafter QSS) are galaxies 
or young globular clusters. Then, a comparatively very 
slow relaxation to thcrmodynamical equilibrium occurs 
— called collisional relaxation — in a timescale of order 
Tcoii ~ N S Td y m where N is the number of particles and 
5 > 1 depends on the system studied. 

The mechanism of collisional relaxation is still not 
well understood. In the context of gravitational systems, 
Chandrasekhar found theoretically, in a seminal work Q , 
an estimate of the relaxation time for gravitational sys- 
tems in three dimensions. He considered an homoge- 
neous system and computed the change in velocity due 
to successive independent collisions of a test particle in 
a stationary macroscopic configuration. Because of the 
hypothesis of homogeneity there is no macroscopic scale 
in the system, which led to an ongoing controversy about 
the value of the maximal impact parameter of the colli- 
sions and in particular how it should scale with N • 
Following this, several studies considered collective ef- 



fects (e.g. 0), but still in homogeneous configurations. 
An explicit theoretical description of the collisional relax- 
ation in inhomogeneous systems is still lacking, despite 
recent progress in this direction 0, Q . 

The collisional relaxation has also been studied numer- 
ically, for a wide variety of systems. In d = 1 dimensions 
for one-dimensional gravity, a scaling of t co u ~ NTd yn has 
been measured for the full relaxation process fioj , and in 
the Hamiltonian Mean Field model the scaling has been 
found to be dependent on the initial condition: t co u 
N 1 - 7 
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[l2|. For dimensions larger than d - 
has been estimated studying — for 
- only its early stage, i.e., for times in which the QSS 
is weakly perturbed (see e.g. [HI 
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, or performing 
simulations with a simplified dynamics. For gravity in 
two-dimensions, in simulations performed imposing ra- 
dial symmetry, it has been observed t co u ~ N 1 - 35 Td yn 
[HI]. In d = 3 dimensions, the Chandrasekhar scaling 
Jroii ~ N/ In N Tdyn has been verified for gravity (e.g. 
0, Ell IHI) and for power-law potential u(r) = 1/r 7 for 
which has been found t co ii ~ Nrdyn if 7 < 1, see [ldllfjj]. 

In this letter, we present an analytical and numeri- 
cal study of the collisional relaxation of a self-gravitating 
system in d = 2 dimensions. The interacting potential 

— solution of the Poisson equation in d — 2 dimensions 

— is u(r) = g\n(r) 7 where g is the coupling constant. It 
is an attractive model to understand the realistic d = 3 
gravity: it presents the same mechanism of collisions as 
in d = 3 (which is not the case for models in d — 1), the 
system is self-confined (it is not necessary to confine it 
artificially in a box), thermal equilibrium properties are 
easily calculated and numerical simulations are easier to 
perform than in d = 3. In this work, we simulated the re- 
laxation dynamics of the system, for the whole time range 
between the QSS and the final thermal equilibrium. The 
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simulations are very well described by a an approximate 
kinetic equation we derive in what follows. 

Theory.- We model the generic evolution of the sys- 
tem using the Boltzmann equation for the one point prob- 
ability density function f(r,v,t). We can write it for- 
mally as 



df 
dt 



df 
Or 



F[f] 



df 
dv 



= r c [/], 



(i) 



where T c [f] is the collision operator. During the relax- 
ation process, the system reaches first a QSS and then 
evolves (comparatively slowly) through an infinity se- 
quence of QSS, in which v • % + F[f] ■ = 0. To 
make Eq. ([T]) tractable analytically, we neglect its sec- 
ond and third term, which implies not taking collective 
effects into account. We focus on the evolution of the ve- 
locity distribution s(v, t) = f d?r /(r, v, t). We integrate 
Eq. (U) over the positions, obtaining, in the approxima- 
tion described above 



ds 

dt 



d 2 rT c [f}. 



(2) 



In the same manner as in the most studied d = 3 case, the 
relaxation is dominated by soft collisions (see e.g. |17j]). 
which implies that the change in velocity of a particle in 
a collision is much smaller than its velocity. Moreover, 
it has been shown that, for times larger than one orbital 
period, the force correlation function decays rapidly (e.g. 
as ~ 1/t 5 for gravity in d = 3 [18]). We may then consider 
that collisions are independent and the use of a Fokker- 
Planck approximation of Eq. (2J is therefore justified (see 
e.g. EHI), which can be written as 



ds(v,t) 
dt 



d_ 

dv,. 
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[D ViVj s(v,t)] 



(3) 

The terms D Vi (v,t) and D ViVj (v,t) are the diffusion co- 
efficients, which are defined as the average change of the 
velocity of the particles per unit of time, (Avi) /At and 
(AviAvj) I At, respectively. 

We estimate first the change in velocity of a particle 
due to a "collision" with a single particle. In the context 
of long range systems, we define a "collision" between two 
particles as the process in which they cross each other in 
one orbital period. The spatial density distribution at 
thermal equilibrium of a self-gravitating system in d = 2 
generates a (mean-field) gravitational potential of the 
form ty(r) — 2 In (A 2 + Ar 2 ) /(3, where /3 is the inverse 
temperature and A and A are two constants which depend 
on the number of particles N and on the total energy of 
the system (see e.g. 15|). For r < A (which corresponds 



to a scale in which are included half of the particles), the 
potential is harmonic, i.e., *ff(r) ~ 4 In A + w 2 r 2 /2 (where 
uj = 4A//3X 2 ). It has been show numerically that this is 
also true in d — 3 for a wide set of initial conditions [2(| ■ 
Assuming that the potential has this approximate form 




FIG. 1: Sketch of the orbits (dotted curves) of two "colliding" 
particles (which are plotted at the same arbitrary time) . The 
plain curve represents their relative trajectory, and the thick 
portion (of length ~ 2b) the part of the trajectory in which 
|AVj_| changes significantly (see text). 



in the QSS (hypothesis which we have checked with our 
numerical simulations), the trajectories of the particles 
in the central region of the system (where collisional re- 
laxation occurs) are ellipses. The relative motion of two 
particles is also therefore an ellipses which can be written 
as r(t) = (xq cos(wi), yo sin(wi)), see Fig. [TJ Because col- 
lisional relaxation is dominated by weak collisions, i.e., 
by the ones in which the trajectory of the particles are 
weakly perturbed, the change in relative velocity in the 
y direction of two crossing particles can be well approxi- 
mated by 



IAV, 



.9 



7r/u 



2/0 dt 



o Xq sin (u)t) + yjj cos 2 (ojt) 



gn 

LUXq 



From geometrical arguments, it is possible to see that 
most of the orbits will have large cllipticity (for example, 
for elliptical orbits with (xo,yo) randomly distributed in 
a finite interval and random phases — which gives rise 
to an approximately homogeneous spatial distribution — 
we find numerically yo/xo ~ 0.3 on average). If we choose 
the axis in order of xq > yo, the distance of closest ap- 
proach is yo and occurs at t = 0. As the integral con- 
verges rapidly an excellent approximation of Q consists 
in taking as upper cutoff of the integral u)t ~ yo/x . If 
yo/xo <C 1, and taking into account that collisions are 
very likely to happen in a central region of the QSS, we 
can therefore approximate the relative velocity during the 
collision by |V(i = 0)| = V ~ lox - Then 



IAV, |~ ^ 



gn 



(4) 



where V is the relative velocity at the distance of closest 
approach. Interestingly, this is the same results than the 
one obtained in the spatial homogeneous case originally 
treated by Chandrasckhar applied to gravity in d = 2. In 
this study, it was considered a rectilinear relative trajec- 
tory with constant relative velocity V (e.g. 17;]), in which 
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the distance of closest approach xq is the impact factor 
b. Then |AVJ ~ 2f™g/[b 2 + (Vt) 2 }dt — which cor- 
responds in this context to the change of velocity in the 
perpendicular direction to the relative velocity — gives 
exactly ^ (see Fig. Q]). This result can be understood 
because an excellent approximation to the latter integral 
is obtained taking t ~ V/b as upper cutoff. Therefore we 
can conclude that the result (HJ) would be a good approx- 
imation when the curvature of the relative trajectories at 
the position of the distance of closest approach is small 
compared with the inverse of this distance. The argu- 
ments presented above apply also in d = 3, which may 
explain why the original Chandrasekhar estimate gives a 
good estimate of the relaxation time in inhomogeneous 
systems, taking as maximal impact parameter the size of 
the system (see e.g. j(3. \\^). 

We see that in the very particular case of two- 
dimensional gravity, the dependence on the impact pa- 
rameter b does not appear in Eq. This implies that, 
in the limit of validity of our approximations, all the par- 
ticles contribute equally to the relaxation, independently 
of their distance. It is possible to compute the change 
in the relative parallel velocity using |AVjJ = VsinO 
and |AV||| = V(l — cos#), where 9 is the angle of de- 
flection. In the weak collision approximation S< 1 and 
thus we have sin.9 ~ 9 and cos6> ~ 1 — 9 2 /2 and then 
|AV||| = |AVj 2 /2y. Taking into account that particle 
masses are equal, we obtain for the change in velocity of 
a particle |Avj_| ^ and |Av||| ~ 

We compute the diffusion coefficients using the stan- 
dard method used in d — 3; we will not detail the calcu- 
lation here. To estimate in our case the number of colli- 
sions per unit of time, we see from our simulations that 
the velocity of the particles are approximately uncorre- 
cted with their position. This allow us not to distinguish 
them. Moreover, we see also that the spatial density dis- 
tribution is approximately constant up to a scale r* (in 
radial coordinates). We can therefore estimate the num- 
ber of collisions of a particle, on average, as NVAt/2r*. 
The mean field (N — > oo) contribution of the change in 
velocity cancels out in average. Averaging over the veloc- 
ity distribution, we obtain, keeping only terms of 0(g 2 ) 
(see [T7| app. L): 



isotropic, we get 



D ViVj (v) 



~aT^ c 

{AvjAvj) 
At 



dgjv) 



(5) 



= C 



d 2 p(v) 
dvidvj 



where C — ir 2 g 2 N/8r* and we have introduced, as in 
the d — 3 case, the Rosenbluth potential [2l[ q(v) 
j d 2 v's(v')/\v — v'\ andp(i>) = J d 2 v's(v')\v — v'| and we 
have assumed that the velocity distribution is isotropic. 
As the succession of QSS have an approximate polar sym- 
metry, it is then useful to write Eq. ((3]) in polar coordi- 
nates. Considering that the Rosenbluth potentials are 



ds 

dt 



c 



d_ 

dv 



q'(v) 



p'(v) 



28V 2 



[p"(v) 



(6) 

We have defined s(v) as the velocity distribution in polar 
coordinates, the primes denotes derivation with respect 
to v and v — |v|/«*. We define the velocity units v* using 
the virial theorem, which states that, for any stationary 
state (and hence a QSS), the average velocity square of 
the particles is constant during all the evolution and it 
is (v 2 ) — gN/2 (e.g. [HI). It is then natural to take 
as velocity unit = y/gN. Equation © depends on N 
through the factor C = C/v%. Therefore C — (Nrdyn)^ 1 , 
which implies t co u ~ Nrdyn, where we have defined the 
dynamical time of the system as Td yn = 1/ ^/gN . 

To compute explicitly the diffusion coefficients we need 
an explicit form of s(v). Taking as velocity distribution 
the equilibrium Gaussian one s(v) — 2vv*f3 exp(— j3v 2 ) 
(with (3 — 2), because of virial theorem, we will make an 
error from only the fourth moment of the velocity. We 
obtain in this approximation 



q(v)=v; 1 e- pv W^PWv 




r,2\ 



(7a) 



e' M +{l + 2Pv 2 )l ((3v 2 ) 



(7b) 



where I n (x) is the modified Bessel function of the first 
kind. It is possible to verify that s(v) = 2vv*f3 exp(— (3v 2 ) 
is a stationary solution of Eq. (|6]) with the diffusion co- 
efficients given by Eq. (J7|). 

Numerical simulations.- We compare the theoret- 
ical model with molecular dynamics simulations per- 
formed with a modification of the publicly available code 
GADGET2 to handle the logarithmic interaction. We 
use a time-step of 2.5 x 10 _4 r c ; yn in order to ensure a 
very precise energy conservation, which is better than 
1CP 5 for the whole duration of the runs. We performed 
simulations with initial water-bag conditions and initial 
virial ratio = ^/2(vq) = 1 and fi = 1.7 (where (vq) 
is the average of the initial velocity square), and differ- 
ent number of particles in the interval N = [100,8000]. 
The simulations have been performed for times between 
t = bOQrdyn for the largest N and t — 750Td yn for the 
smallest ones. In order to improve statistics, we average 
the measured velocity distribution over 100 consecutive 
snapshots in an interval of 2.5Td yn - We monitor how the 
system approaches thermal equilibrium using the param- 
eter £(t) = ^ f™[N(v,t)- N MB (v)] 2 dv, where N MB (v) 
is the thermal equilibrium velocity distribution. In order 
to compare simulations with theory we compute the as- 
sociated Langevin equation of Eq. Therefore, the 
change in the velocity is given, following the Ito defini- 
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FIG. 2: Upper curves: initial condition with y, = 1. Lower 
curves: initial conditions with /j, — 1.7. Points: evolution 
of the crossover parameter £(t) measured in the molecular 
dynamics simulations for the two different initial condition 
fi — 1 and n — 1.7. Lines: theoretical prediction calculated 
using Eq. ((U for each case (see text). 

tion, by 

dv {t) = c\^q\v) + 1 ^^dt+^WW)dwY (8) 

where dW is a Gaussian stochastic variable delta corre- 
lated in time with variance unity. We choose as initial 
condition a configuration of the numerical simulation at 
t/Tdyn = 260 (time in which the system has violently re- 
laxed) and then we compare the evolution predicted by 
the Langevin equation and the one of the full numerical 
simulation. We integrate Eq. (|8]) by a simple Euler pro- 
cedure. In Fig.[2]we show the evolution of £(t), where the 
time axis has been rescaled by a factor N, which indicates 
a scaling of the relaxation time as t co u ~ NTd yn - For 
clarity, between all the simulations with different num- 
bers of particles performed we plot three of them. The 
part of the curve which flattens corresponds to thermal 
equilibrium, which is attained first as N decreases. The 
matching between the curves corresponding to different 
N is very good in the region out of equilibrium, as it has 
been illustrated for N = 750, N = 12 3 and N = 16 3 , 
which confirm the prediction of Eq. ([6]) for the scaling 
of the relaxation. The dashed curves corresponds to the 
theoretical prediction given by Eq. (jSJ with r* — 0.38 for 
the simulation with ji = 1 and r* = 0.2 for the simula- 
tion with fx = 1.7. These values are close to the scale of 
the falloff in the density distribution; the density decays 
to half its center value around r ps 0.4 for both set of 
simulations. Note that the difference in the slopes of the 
curves is essentially due to the different initial conditions 
considered for each case rather than in the value of r* 



FIG. 3: Evolution of the velocity distribution function (left 
panel initial conditions fi = 1 and right panel initial conditions 
fi = 1.7) for the molecular dynamics simulation (straight red 
line) and the prediction of the Langevin equation (dotted blue 
line). The time increases as the shift in the y axis increases. 
At the last time the system has reached thermal equilibrium. 
The pink dashed line (plotted for the first time) corresponds 
to the distribution of thermal equilibrium. 



taken. The full simulation curves decay to a lower value 
at thermal equilibrium because fluctuations appears to 
be larger in the molecular dynamics simulations than in 
the Langevin simulation. In Fig.[3]we show the evolution 
of the full velocity distribution function for both the sim- 
ulation and the theory. We see a very good agreement 
between them, and specially in the [i = 1.7 simulation, in 
which the Langevin equation is able to reproduce prop- 
erly the relaxation of the core-halo structure in velocity 
space. 

Conclusions. - We do not observe the scaling t co u ~ 
^V 1 ' 35 Tdi/n observed in [IBj]. This is is due to the fact that 
their simplified dynamics appear not to describe properly 
the collisional dynamics of the real d — 2 system. 

Some final remarks can be made about the maximal 
impact parameter which has to be considered in the cal- 
culations. In our calculations we don't use any artificial 
spatial cutoff in counting the number of collisions per unit 
of time. This is equivalent to not introducing any cut- 
off in the maximal impact parameter for interactions in 
which the change in velocity depends on it — such as the 
gravity in d — 3. We obtained a good agreement between 
theory and simulations not only in the scaling of t co u 
with N but also in the amplitude of which allows 
us to conclude that it is not necessary to introduce any 
artificial cutoff in the calculations but integrals are auto- 
matically regularized by the size of the system. This re- 
sult is in agreement with simulations performed in d = 3 
dimensions [16J with potential interactions u(r) ~ 1/V 7 , 
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7 < 1, in which the maximal impact parameter to take 
in the Chandrasekhar approximation was numerically es- 
timated to be 1/3 the size of the system. 
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